Render this, for instance, like:
templater::render_template("summarize_run.Rmd",output="lostruct_results/type_snp_size_10000_jobid_324902/run_summary.html",change.rootdir=TRUE)
This run had these parameters:
Here are the number of windows per chromsome, and the computed MDS coordinates, colored by chromosome:
table(regions$chrom)
##
## chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8
## 497 433 542 541 449 365 476 437
pairs( mds.coords[,mds.cols], pch=20, col=adjustcolor(chrom.cols,0.75) )
if (do.pdfs) { pdf_copy() }
## pdf version at: figure/run-summary/mds_pairplot-1.pdf
## pdf
## 2
Here are the extreme windows in the MDS plot:
mds.corners <- corners( mds.coords[,mds.cols[1:2]], prop=.05 )
# set up colors and pchs for corners
corner.cols <- brewer.pal(3,"Dark2")
corner.pch <- c(15,17,19)
ccols <- rep("black",nrow(mds.coords))
cpch <- rep(20,nrow(mds.coords))
for (k in 1:ncol(mds.corners)) {
ccols[ mds.corners[,k] ] <- corner.cols[k]
cpch[ mds.corners[,k] ] <- corner.pch[k]
}
if (do.pdfs) { pdf_copy() }
## pdf version at: figure/run-summary/get_corners-1.pdf
## pdf
## 2
# plot corners and MDS along the chromosome
spacing <- 1
opar <- par(mar=c(4,4,2,1)+.1,mgp=c(2.5,0.8,0))
layout(matrix(c(rep(1,length(mds.cols)),1+seq_along(mds.cols)),ncol=2),
widths=c(1,2), heights=layout_heights(length(mds.cols),dl=spacing,ncol=2))
plot( mds.coords[,mds.cols[1:2]], pch=cpch,
col=adjustcolor(ccols,0.75), asp=1,
xlab="MDS coordinate 1", ylab="MDS coordinate 2" )
opar2 <- par(mar=c(par("mar"),spacing/2)[c(5,2,3,4)])
for (k in mds.cols) {
lastone <- (k==mds.cols[length(mds.cols)])
if (lastone) { par(mar=c(par("mar"),opar2$mar[1])[c(5,2,3,4)]) }
chrom.plot( mds.coords[,k], pch=20,
xlab=if (lastone) { "Position (Mb)"} else { "" }, # main=paste("MDS coordinate",match(k,mds.cols)),
chrom.labels=lastone,
ylab=colnames(mds.coords)[k],
col=adjustcolor(ccols,0.75) )
# do this for all but first
par(mar=c(par("mar"),spacing/2)[c(1,2,5,4)])
}
par(opar)
if (do.pdfs) { pdf_copy() }
## pdf version at: figure/run-summary/plot_corners-1.pdf
## pdf
## 2
Now, we'll look at PCA plots from the extracted corners. (this is done without a ton of memory by accumulating the covariance matrix in running_cov):
corner.npc <- 4
corner.regions <- lapply( 1:ncol(mds.corners), function (k) {
regions[ mds.corners[,k],]
} )
corner.covmats <- lapply( 1:ncol(mds.corners), function (k) {
reg <- regions[ mds.corners[,k], ]
qfun <- multi_vcf_query_fn( chrom.list=chroms, file=bcf.files, regions=reg )
running_cov(qfun,1:nrow(reg))
} )
corner.pca <- lapply( corner.covmats, function (covmat) {
cov_pca(covmat=covmat,k=corner.npc,w=opt$weights)
} )
Here is the color scheme:
pop.names <- levels(samps$population)
pop.cols <- rainbow_hcl(nlevels(samps$population))
pop.pch <- seq_len(nlevels(samps$population))
plot( rep(1,length(pop.names)), seq_along(pop.names), pch=pop.pch, col=pop.cols, xlim=c(0,length(pop.names)),
xlab='', ylab='', xaxt='n', yaxt='n' )
text( rep(1,length(pop.names)), seq_along(pop.names), labels=pop.names, pos=4 )
Here are all pairwise plots of the first 4 PCs for each of the three corners:
layout(t(1:3))
for (i in 1:(corner.npc-1)) {
for (j in (i+1):corner.npc) {
for (k in 1:ncol(mds.corners)) {
vectors <- matrix( corner.pca[[k]][-(1:(1+corner.npc))], ncol=corner.npc )[,c(i,j)]
colnames(vectors) <- paste("PC", c(i,j))
plot(vectors, pch=pop.pch[samps$population],
col=pop.cols[samps$population] )
if (i==1 && j==2) {
mtext(paste("corner",k),side=3,cex=2)
}
}
if (do.pdfs) { pdf_copy(plot.id=paste(i,j,sep="_")) }
}
}
## pdf version at: figure/run-summary/plot_corner_pca-1_1_2.pdf
## pdf version at: figure/run-summary/plot_corner_pca-1_1_3.pdf
## pdf version at: figure/run-summary/plot_corner_pca-1_1_4.pdf
## pdf version at: figure/run-summary/plot_corner_pca-1_2_3.pdf
## pdf version at: figure/run-summary/plot_corner_pca-1_2_4.pdf
## pdf version at: figure/run-summary/plot_corner_pca-1_3_4.pdf